function solve = solve(j,fp,param)
alpha = param(1:9);
TFP = param([10:11 end]) ;  


h =1; %solve for the equilibrium first;
College = (fp.College_forward{j}(:,h)<=fp.schools_distrib(j,4));

fp_parfor = struct();
fp_parfor.mean_skills    = fp.schools_distrib(j,1).*(College==0) + fp.schools_distrib(j+4,1).*(College==1);  % # mean skills (by mother education); 
fp_parfor.sd_skills      = fp.schools_distrib(j,2).*(College==0) + fp.schools_distrib(j+4,2).*(College==1);  % # standard deviation skills (by mother education);
fp_parfor.n_simulation   = fp.schools_distrib(j,3);  % # children in the economy;


fp_parfor.init_draws = fp_parfor.mean_skills + fp_parfor.sd_skills.*(fp.init_draws{j}(1:fp_parfor.n_simulation,1,h) - mean(fp.init_draws{j}(1:fp_parfor.n_simulation,1,h)))./std(fp.init_draws{j}(1:fp_parfor.n_simulation,1,h));    
fp_parfor.peers_forward  = fp.peers_forward{j}(1:fp_parfor.n_simulation,1:fp_parfor.n_simulation,h,:);
fp_parfor.peers_backward = fp.peers_backward{j}(:,1:fp_parfor.n_simulation,:,:);    
fp_parfor.PS_forward     = fp.PS_forward{j}(1:fp_parfor.n_simulation,h,:);
fp_parfor.College    = College;

epsilon = 5;
toll= 0.05;

epsilon_skills = 5;

record_epsilon=[3 epsilon];

i = 1; %iteration of the fixed point;

while (epsilon > toll  && epsilon_skills>toll) %Fix point on the policy functions;
    
fp_parfor.i = i; %iteration   
fp_parfor.j = j; %neighborhood      
fp_parfor.h = h; %simulation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;    
%Simulate Forward;    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

if i == 1

    %first iteration: set the initial guess about the distribution of
    %skills;
    fp_parfor.distrib_PS = zeros( size( fp_parfor.init_draws , 1 ) ,fp.periods-1);     
    fp_parfor.inv0{1} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    fp_parfor.inv0{2} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    fp_parfor.inv0{3} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    
    fp_parfor.inv1{1} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    fp_parfor.inv1{2} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    fp_parfor.inv1{3} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);    
    
    distrib_skills1 = exp(fp_parfor.init_draws(:,1));            
    distrib_skills2 = technology(1,alpha, distrib_skills1 , fp.Max_Time, fp_parfor.distrib_PS(:,1) ,  distrib_skills1 , fp_parfor.College, TFP) ;
    distrib_skills3 = technology(2,alpha, distrib_skills2 , fp.Max_Time, fp_parfor.distrib_PS(:,2) ,  distrib_skills2 , fp_parfor.College, TFP)  ;
    distrib_skills4 = technology(3,alpha, distrib_skills3 , fp.Max_Time, fp_parfor.distrib_PS(:,3) ,  distrib_skills3 , fp_parfor.College, TFP)  ;        
    
    fp_parfor.distrib_skills(:,1) = distrib_skills1;
    fp_parfor.distrib_skills(:,2) = distrib_skills2;
    fp_parfor.distrib_skills(:,3) = distrib_skills3;
    fp_parfor.distrib_skills(:,4) = distrib_skills4;
    

    fp_parfor.grid{1} = prctile( distrib_skills1 , fp.percentiles_kid ) ;
    fp_parfor.grid_peers{1} = fp_parfor.grid{1} ;    
    fp_parfor.grid{2} = prctile( distrib_skills2 , fp.percentiles_kid ) ;
    fp_parfor.grid_peers{2} = fp_parfor.grid{2} ;              
    fp_parfor.grid{3} =  prctile(distrib_skills3 , fp.percentiles_kid );
    fp_parfor.grid_peers{3} = fp_parfor.grid{3};

else
    %solve forward to simulate the distribution of skills and investments;
    sf = solve_forward(param,fp,fp_parfor,h);
        
    previous_iter_distrib_skills = fp_parfor.distrib_skills;    
    
    fp_parfor.distrib_skills = sf.distrib_skills;
    fp_parfor.distrib_PS = sf.distrib_PS;
    
    fp_parfor.grid{1} = linspace( prctile(exp(fp_parfor.init_draws(:,1)) , 5) , prctile(exp(fp_parfor.init_draws(:,1)), 95), length(fp.percentiles_kid) );
    fp_parfor.grid_peers{1} = fp_parfor.grid{1} ;
    fp_parfor.grid{2} = linspace( prctile( fp_parfor.distrib_skills(:,2) , 5) , prctile(fp_parfor.distrib_skills(:,2), 95 ), length(fp.percentiles_kid) );
    fp_parfor.grid_peers{2} = fp_parfor.grid{2} ;              
    fp_parfor.grid{3} = linspace( prctile( fp_parfor.distrib_skills(:,3), 5 ) , prctile(fp_parfor.distrib_skills(:,3), 95 ), length(fp.percentiles_kid) );
    fp_parfor.grid_peers{3} = fp_parfor.grid{3} ;    
 
   fp_parfor.inv0{1} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
   fp_parfor.inv0{2} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
   fp_parfor.inv0{3} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    
   fp_parfor.inv1{1} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
   fp_parfor.inv1{2} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
   fp_parfor.inv1{3} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);      

    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Solve Backward;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

sb= solve_backward(param, fp,fp_parfor) ;      

fp_parfor.Hc_tp1_guess_prime = sb.Hc_tp1_guess_prime;

if i > 1 
eps1=zeros(fp.periods-1,1);
eps2=zeros(fp.periods-2,1);
eps_skills=zeros(fp.periods-2,1);


for t = 1:fp.periods-1
            
  eps1(t) = max(abs(fp_parfor.estim_policy0{t}(:) - sb.estim_policy0{t}(:)));
  
  if t==3

  end
  
    if t < fp.periods-1
        
        eps2(t) = max(abs(fp_parfor.estim_policy1{t}(:) - sb.estim_policy1{t}(:)));
 
    end
    
   if t > 1  
   deciles_skills_previous_iter = quantile(log(previous_iter_distrib_skills(:,t)),[0.05:0.1:0.95]);
   deciles_skills_current_iter = quantile(log(sf.distrib_skills(:,t)),[0.05:0.1:0.95]);          
   
   eps_skills(t) =  max(abs(deciles_skills_previous_iter-deciles_skills_current_iter));
          
   end    
    
    
end

epsilon = max([eps1;eps2]);
epsilon_skills = max(eps_skills);
%keep track of epsilon records;
ind_rec = mod(i-1,2) + 1;
record_epsilon(ind_rec) = epsilon;


end




fp_parfor.estim_policy0 = sb.estim_policy0;
fp_parfor.estim_policy1 = sb.estim_policy1;

fp_parfor.estim_value0 = sb.estim_value0;
fp_parfor.estim_value1 = sb.estim_value1;

if i ==1
fprintf('Solving Neighborhood %d:', j)
end
fprintf('%d ', i)

i = i + 1;
%record_epsilon


end

fprintf('Solved Neighborhood %d \n ', j)

Value0(:,:,:,:,:,h) = sb.Value0_actual;
Value1(:,:,:,:,:,h) = sb.Value1_actual;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Simulate forward once the fixed point is achieved and create the simulated 
%dataset;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

fprintf('Simulating Neighborhood %d:\n ', j )

for h = 1 : fp.n_rip_sim %Loop over several simulation of the economy;

if h ==fp.n_rip_sim %Last simulation for this neighborhood;
fprintf('Simulated Neighborhood %d \n ', j)
end

College = (fp.College_forward{j}(:,h)<=fp.schools_distrib(j,4));
fp_parfor.mean_skills    = fp.schools_distrib(j,1).*(College==0) + fp.schools_distrib(j+4,1).*(College==1);  % # mean skills (by mother education); 
fp_parfor.sd_skills      = fp.schools_distrib(j,2).*(College==0) + fp.schools_distrib(j+4,2).*(College==1);  % # standard deviation skills (by mother education);

fp_parfor.init_draws     = fp_parfor.mean_skills + fp_parfor.sd_skills.*(fp.init_draws{j}(1:fp_parfor.n_simulation,1,h) - mean(fp.init_draws{j}(1:fp_parfor.n_simulation,1,h)))./std(fp.init_draws{j}(1:fp_parfor.n_simulation,1,h));
fp_parfor.College    = College;

fp_parfor.peers_forward  = fp.peers_forward{j}(1:fp_parfor.n_simulation,1:fp_parfor.n_simulation,h,:);
fp_parfor.PS_forward     = fp.PS_forward{j}(1:fp_parfor.n_simulation,h,:);    
    
fp_parfor.i = i;    
fp_parfor.j = j;
fp_parfor.h = h;

%Generate simulated data;
sf = solve_forward(param,fp,fp_parfor,h);

%Store data for each simulation;
sim_data(:,:,h) = sf.sim_data;
sim_data_network(:,:,h) = sf.sim_data_network;

end

solve.Value0 = Value0;
solve.Value1 = Value1;
solve.sim_data = sim_data;
solve.sim_data_network = sim_data_network;